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We consider the trajectories of particles diffusing between two infinite baths of fixed concentrations 
connected by a channel, e.g. a protein channel of a biological membrane. The steady state influx 
and efflux of Langevin trajectories at the boundaries of a flnite volume containing the channel and 
parts of the two baths is replicated by termination of outgoing trajectories and injection according 
to a residual phase space density. We present a simulation scheme that maintains averaged flxed 
concentrations without creating spurious boundary layers, consistent with the assumed physics. 

PACS numbers: 83.10.Mj, 02.50.-r, 05.40.-a. 



I. INTRODUCTION 

We consider particles that diffuse in a domain fl 
connecting two regions, where fixed, but possibly 
different concentrations are maintained by connec- 
tion to practically infinite reservoirs. This is the 
situation in the diffusion of ions through a protein 
channel of a biological membrane that separates 
two salt solutions of different fixed concentrations 

El- 

Continuum theories of such diffusive systems de- 
scribe the concentration field by the Nernst-Planck 
equation (NPE) with fixed boundary concentra- 
tions 01-01 • On the other hand, the underlying 
microscopic theory of diffusion describes the mo- 
tion of the diffusing particles by Langevin's equa- 
tions i3, This means that on a microscopic 
scale there are fluctuations in the concentrations 
at the boundaries. The question of the boundary 
behavior of the Langevin trajectories (LT), corre- 
sponding to fixed boundary concentrations, arises 
both in theory and in the practice of particle sim- 
ulations of diffusive motion 

When the concentrations are maintained by con- 
nection to infinite reservoirs, there are no physical 
sources and absorbers of trajectories at any definite 
location in the reservoir or in fl. The boundaries 
in this setup can be chosen anywhere in the reser- 
voirs, where the average concentrations are effec- 
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tively fixed. Nothing unusual happens to the LT 
there. Upon reaching the boundary they simply 
cross into the reservoir and may cross the bound- 
ary back and forth any number of times. Limit- 
ing the system to a finite region necessarily puts 
sources and absorbers at the interfaces with the 
baths, as described in flE\. 

The boundary behavior of diffusing particles in a 
finite domain H. has been studied in various cases, 
including absorbing, refiecting, sticky boundaries, 
and many other modes of boundary behavior [T^ , 
tlZj . fn a sequence of Markovian jump pro- 
cesses is constructed such that their transition 
probability densities converge to the solution of the 
Nernst-Planck equation with given boundary con- 
ditions, including fixed concentrations and sticky 
boundaries. Brownian dynamics simulations with 
different boundary protocols seem to indicate that 
density fluctuations near the channels are indepen- 
dent of the boundary conditions, if the boundaries 
are moved sufficiently far away from the channel 
iil9>i . However, as shown in |,2Q.j , many bound- 
ary protocols for maintaining flxed concentrations 
lead to the formation of spurious boundary layers, 
which in the case of charged particles may produce 
large long range fluctuations in the electric held 
that spread throughout the entire simulation vol- 
ume ri. The analytic structure of these boundary 
layers was determined in |2^. following several 
numerical investigations (e.g, [23j)- 

ft seems that the boundary behavior of LT of 
particles diffusing between flxed concentrations has 
not been described mathematically in an adequate 
way. From the theoretical point of view, the ab- 
sence of a rigorous mathematical theory of the 
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boundary behavior of LT diffusing between fixed 
concentrations, based on the physical theory of the 
Brownian motion, is a serious lacuna in classical 
physics. 

It is the purpose of this letter to analyze the 
boundary behavior of LT between fixed concen- 
trations and to design a Langevin simulation that 
does not form spurious boundary layers. We find 
the joint probability density function of the veloci- 
ties and locations, where new simulated LT are in- 
jected into a given simulation volume, while main- 
taining the fixed concentrations. As the time step 
decreases the simulated density converges to the 
solution of the Fokker-Planck equation (FPE) with 
the imposed boundary conditions without forming 
boundary layers. 



II. TRAJECTORIES, FLUXES, AND 
BOUNDARY CONCENTRATIONS 

We assume fixed concentrations Cl and Cji on 
the left and right interfaces between Vl and the 
baths i?, respectively, with all other boundaries of 
VL being impermeable walls, where the normal par- 
ticle flux vanishes. We assume that the particles 
interact only with a mean field, whose potential is 
$(a;), so the diffusive motion of a particle in the 
channel and in the reservoirs is described by the 
Langevin equation (LE) 



x{0) = xo, v{0) = vo. 



(1) 



where 7(3;) is the (state-dependent) friction per 
unit mass, e is a thermal factor, and ti; is a vector of 
standard independent Gaussian (5-correlated white 
noises ^. 

The probability density function (pdf) of finding 
the trajectory of the diffusing particle at location x 
with velocity v at time t, given its initial position, 
satisfies the Fokker-Planck equation (FPE) in the 
bath and in the reservoirs. 



dp 



= -v-\7xP + lix)eAvP (2) 



-V 



V ■ 



7(a;)t> + Va;$(a;) 



p{x,v,0\xo,vo) = S{x - Xo,V - Vo). 

In the Smoluchowski limit of large friction the sta- 



tionary solution of (|2Jl admits the form |^ 
where the flux density vector J(x) is given by 



(3) 



J{x) = --l-^S^eWp{x)+p{x)W<i>{x)'^+0 (^-^ 



and p{x) satisfies 



-V • J{x) = V 



7 (a;) 



eVp{x) +p{x) V$ (x) 



In one dimension, the stationary pdfs of velocities 
of the particles crossing the interface into the given 
volume are 



Vl{v) 



Jv 




for 7; > 0, 



(4) 



for w < 0, 



where J is the net probability flux through the 
channel. The source strengths (unidirectional 
fluxes at the interfaces) are given by |^ 



Jl = xI—Cl-^ + 0{- 
J 



(5) 



Ce 
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III. APPLICATION TO SIMULATION 

Langevin simulations of ion permeation in a pro- 
tein channel of a biological membrane have to 
include a part of the surrounding bath, because 
boundary conditions at the ends of the channel 
are unknown. The boundary of the simulation has 
to be interfaced with the bath in a manner that 
does not distort the physics. This means that new 
LT have to be injected into the simulation at the 
correct rate and with the correct distribution of 
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displacement and velocity, for otherwise, spurious 
boundary layers will form poj . 

Consider a single simulated trajectory that 
jumps according to the discretized LE ^ 

x{t + At) = v{t)At, v{t + At) = v{t){l - -/At) 
- Vx^ixit))At + y^Aw{t), (6) 

where Aw is normally distributed with zero mean 
and variance At. The trajectory is terminated 
when it exits f2 for the first time. The problem 
at hand is to determine an injection scheme of 
new trajectories into f2 such that the interface con- 
centrations are maintained on the average at their 
nominal values C'l and Cr and the simulated den- 
sity profile satisfies 

To be consistent with |(2Jl, the injection rate has 
to be equal to the unidirectional flux at the bound- 
ary (|SJ). New trajectories have to be injected with 
displacement and velocity as though the simulation 
extends outside fl, consistently with the scheme 
©, because the interface is a fictitious boundary. 
The scheme ^ can move the trajectory from the 
bath B into f2 from any point ^ € B and with 
any velocity rj. The probability that a trajectory, 
which is moved with time step At from the bath 
into Q, or from into the bath will land exactly 
on the boundary is zero. It follows that the pdf of 
the point (a;, v), where the trajectory lands in Q in 
one time step, at time t' = t + At, say, given that 
it started at a bath point (^, rj) (in phase space) 
is, according to ©, 



Pr{x{t') ^ X, v{t') = V I x{t) = v{t) = ?7} 
5{x-i- r]At) 



exp ■ 



(47r£7At)3/2 

\v-r]- {jv + V<^{^))At\ 
4ejAt 



(7) 



-o(At). 



Thus the first point of a new trajectory is chosen 
according to the pdf ^ and the subsequent points 
are generated according ((HJ, that is, with the tran- 
sition pdf lO, until the trajectory leaves fl. By 
construction, this scheme recovers the joint pdf (|2Jl 
in r2, so no spurious boundary layer is formed. 

As an example, we consider a one-dimensional 
Langevin dynamics simulation of diffusion of free 
particles between fixed concentrations on a given 
interval. Assuming that in a channel of length L 



The stationary pdf p{^, rj) of such a bath point is 
given in Q . The conditional probability of such a 
landing is 

Fr{x,v\xen,^e B}^ (8) 

^3 dr, ^ d| Pr{vit') - V, xit') =x\t v}pit v) 

Fr{x £n,^eB} ' 

where the denominator is a normalization constant 
such that 

^dv dxPr{x,v\x en,$,eB} = 1. 



which means that 7 is sufficiently large, the flux 
term in eq.® is negligible relative to the concen- 
tration term. The concentration term is linear with 
slope and thus can be approximated by a con- 
stant, so that p{i) = p(0) -I- O (7-1) in the left 
bath. Actually, the value of p(0) 7^ is unimpor- 
tant, because it cancels out in the normalized pdf 
IHl), which comes out to be 



exp 



¥v{x,v\x > 0,$ < 0} = 



2e[l + {^Atf 



2eAt^l + (7At)2 



(9) 



4e7At \At ^ + (7Ai)2 



In the limit At ^ Q we. obtain from cq.Q 



¥v{x,v\x > 0,^ < 0} 



V27re 



where H(v) is the Heaviside unit step function. 
This means that with the said approximation, LT 
enter at a; = with a Maxwellian distribution of 
positive velocities. Without the approximation the 
limiting distribution of velocities is Note, how- 
ever, that injecting trajectories by any Markovian 
scheme, with the limiting distribution and 
with any time step At, creates a boundary layer 

A LD simulation with Cl 7^ 0, Cr = 0, and the 
parameters 7 — 100, e — 1, L ~ 1, At = 10~^ 
with 25000 trajectories, once with a Maxwellian 
distribution of velocities at the boundary x = 
(red) and once with the pdf © (blue) shows that 
a boundary layer is formed in the former, but not 
in the latter (see Figure 

An alternative way to interpret eq. Q is to view 
the simulation jnj as a discrete time Markovian 
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FIG. 1: Left panel: Concentration against displace- 
ment of a LD simulation with injecting particles ac- 
cording to the residual distribution ^ (blue), and ac- 
cording to the Majcwellian velocity distribution IIUI 
exactly at the boundary (red). The two graphs are al- 
most identical, except for a small boundary layer near 
a; = in red. Right panel: Zoom in of the concentra- 
tion profile in the boundary layer x < 0.01 = •s/e/7. 



process {x{t),v{t)) that never enters or exits fl 
exactly at the boundary. If, however, we run a 
simulation in which particles are inserted at the 
boundary, the time of insertion has to be random, 
rather than a lattice time nAt. Thus the time of 



the first jump from the boundary into the domain 
is the residual time At' between the moment of in- 
sertion and the next lattice time {n + l)At. The 
probability density of jump size in both variables 
has to be randomized with At', with the result @. 
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